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Abstract: A parametric level set method (PaLS) is implemented for 
image reconstruction for hyperspectral diffuse optical tomography (DOT). 
Chromophore concentrations and diffusion amplitude are recovered using 
a linearized Born approximation model and employing data from over 100 
wavelengths. The images to be recovered are taken to be piecewise constant 
and a newly introduced, shape-based model is used as the foundation 
for reconstruction. The PaLS method significantly reduces the number of 
unknowns relative to more traditional level-set reconstruction methods and 
has been show to be particularly well suited for ill-posed inverse problems 
such as the one of interest here. We report on reconstructions for multiple 
chromophores from simulated and experimental data where the PaLS 
method provides a more accurate estimation of chromophore concentrations 
compared to a pixel-based method. 
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1. Introduction 

Near-infrared light has proven to be useful for imaging the human body and providing func- 
tional information for applications including breast cancer detection and characterization [1-5]. 
In recent years significant improvements have been made in diffuse optical tomography (DOT) 
by using data from multiple wavelengths. The use of this kind of hyperspectral information 
has moved DOT away from the recovery of space and time varying maps of absorption and 
scattering properties of the breast to the recovery of chromophore concentrations as well as 
diffusion amplitude. Inverting for multiple chromophores such as hemoglobin, lipids and water 
has been shown to provide an improved ability to localize tumours or other objects of interest 
in the breast cancer application [6-8]. While these and other advancements have been critical 
for moving DOT from the lab into the clinic, there remain significant obstacles to be addressed 
so that the method can be used to stably recover these quantities. 

The imaging of chromophore concentration and diffusion amplitude from DOT data requires 
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the solution of an ill-posed and non-linear inverse problem [3]. The physics associated with 
the photons travelling through turbid media, such as breast tissue, makes the problem ill-posed 
and therefore introduces significant challenges for the DOT problem. The ill-posedness make 
the reconstruction highly sensitive to noise and un-modeled effects which can reduce accuracy 
in the recovered images [7]. Additional problems are encountered when recovering multiple 
chromophores such as crosstalk, e.g. when two spatially disjoint chromophores pollute the 
images of each other. 

The ill-posedness of the DOT problem is traditionally solved by putting the recovery of 
images in a variational context. To achieve this, numerous regularization techniques [9] can 
be used to reduce artifacts and improve the overall accuracy. Traditionally, methods such as 
Tikhonov and total variation (TV) regularization have been employed along with L-curve meth- 
ods to optimally choose regularization parameters [10-12]. Increasing the amount of data used 
to solve the DOT problem has also proven effective in generating high accuracy images [8, 13]. 
To this end hyperspectral information is employed, which along with regularization has been 
shown to reduce the non-uniqueness of the solution to the image formation problem [8, 10, 1 1]. 

In this paper, we expand on our previous work in [10] where we employed hyperspectral 
information for recovery of chromophore concentrations. Implementing the Born approxima- 
tion along with Tikhonov regularization we demonstrated that hyperspectral information was 
shown to have utility for a pixel -based reconstruction of chromophores [10]. The large data set 
introduced significant challenges in constructing and computing with the forward model and 
great care had to be taken in choosing the regularization parameters. The known limitations 
of the Born approximation and the high number of unknowns encountered in the pixel-based 
approach, made reconstructions sensitive to the choice of regularization parameters [14]. This 
motivated us to move to a geometric reconstruction approach based on a low-order parametric 
model that has been shown to perform well in the context of ill posed inverse problems [15-17]. 

The primary contribution of this work is the development of a new approach to the multi- 
chromophore inverse problem in which the medium is or can be well approximated as piecewise 
constant. Piecewise constant DOT problems have been considered in the past. In [18] level sets 
were used in a two-step method for shape estimation assuming that prior information of the 
absorption parameter was known. Schweiger et al. [19] and Kilmer et al. [20] employed level 
sets for the DOT problem estimating parameter distributions using a piecewise basis. Arridge 
et al. investigated shape based methods by estimating level-sets, specifically investigating an 
explicit method using basis functions and an implicit shape reconstruction to recover absorption 
and diffusion assuming a known background [21]. For a detailed review of the use of level sets 
in inverse scattering problems we point the reader to [22]. 

The approach we consider in this paper is significantly different from those in [18, 19, 21]. 
In addition to the fact that none of these papers have considered the fully hyperspectral case, 
some of these methods [18-20] require the recovery of unknown quantities defined on a fine 
scale pixelated discretization of the region of interest. More specifically in [21] absorption and 
scattering is estimated using level sets assuming the background known. With the Born ap- 
proximation we assume the absorption and diffusion is known in the background but here we 
estimate chromophore concentrations and diffusion of the object of interest as well as chro- 
mophore concentration in the background. TV reconstructions recover images while traditional 
level set methods work with a level set function defined on a pixel-based grid. In both cases, 
regularization is required to obtain adequate results and one is faced with the corresponding 
challenge of choosing regularization parameters [1 1, 19]. 

In this paper we consider the use of a shape-based approach to the hyperspectral DOT prob- 
lem based on a newly-developed parametric level set (PaLS) formulation. In [15], a basis func- 
tion expansion was used to provide a low order representation of the level set function and 



#162145 - $15.00 USD 
(C) 2012 OSA 



Received 27 Jan 2012; accepted 15 Mar 2012; published 18 Apr 2012 
1 May 2012 / Vol. 3, No. 5 / BIOMEDICAL OPTICS EXPRESS 1009 



yielded very strong results for a number of highly ill-posed inverse problems including a re- 
stricted form of the DOT problem where a single wavelength was employed to determine only 
optical absorption. The method in [15] required no explicit regularization and, due to the low- 
order nature of the model (number of parameters much, much less than number of pixels) was 
amenable to Newton-type inversion algorithms known to converge more rapidly than gradient- 
based schemes. Moreover, in [15] our experiments indicated a surprising lack of sensitivity to 
the initialization of the inversion algorithm. 

Although the breast is a highly heterogeneous medium, in the level set method we assume the 
images to be recovered to be piecewise constant. This approximation is supported in the litera- 
ture. For example Schweiger et al. assumed anatomical prior information to derive a piecewise 
constant region basis [4]. The co-located structures of different species of chromophores is dis- 
cussed further in Section 5. In the future, where we move to more complicated simulations and 
experiments the goal would be to implement multiple characteristic functions that require the 
estimation of multiple level sets, discussed further in Section 8. To prove the accuracy of our 
method, we use simulation data to show parametric level-set reconstructions for multiple chro- 
mophores and diffusion amplitude. Experimental data are used to show the improvement of the 
PaLS method over pixel based methods using previously reported data sets [10]. 

The remainder of this paper is structured as follows. In Section 2 we discuss our forward 
model and derive in detail the Born model used. In Section 3 we discuss the PaLS method and 
in Section 4 we derive the Levenberg Marquardt algorithm used for the inverse problem. In 
Section 7.1 we present simulation results with multiple chromophore and diffusion amplitude 
reconstructions, and in Section 7.2 we present experimental results that show the improvement 
of the PaLS method over pixel based methods. 

2. Forward problem 

In this paper we restrict our attention to problems in which the transport physics [2] associated 
with the propagation of light at wavelength A through tissue can be adequately approximated 
using a diffusion model [23,24] of the form 

V ■D 0 (r,X)V<i>(r,X)+vrf(r,X)<P(r,X) = -vS(r,X) (1) 

where <J>(r,A) is the photon fluence rate at position r due to light of wavelength A injected 
into the medium, v is the electromagnetic propagation velocity in the medium, n®(r,X) is the 
spatially varying absorption coefficient, and S(r,X) is the photon source with units of optical 
energy per unit time per unit volume. For the work in this paper the sources are considered to 
be delta sources in space and can be written as S(r, X) = So(X)8(r — r s ) with 5b(A) the source 
power at wavelength X. For spatially varying scattering we assume that the diffusion coefficient 
D°(r,X) follows Mie scattering theory where a scattering prefactor *P depends on the size and 
density of scatterers while a scattering exponent b depends on the size of the scatterers. Using 
this, we write the perturbation in the diffusion coefficient as 

^) = 3^(£)* = VA ^W- (2) 

The arbitrarily chosen reference wavelength Ao is introduced to achieve a form of the Mie 
model where *P has the units of mm -1 and X P' has units of mm. In this paper, for simplicity we 
consider the case where the background medium is infinite and homogeneous. Generalization 
to the more practical case where there are boundaries is straightforward in theory though some- 
what more complex in terms of implementation [25]. As the primary objective of this work is 
the demonstration of a new approach to inversion, we prefer to consider the simpler physical 
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problem holding off on the more realistic implementation for the future. Additionally, we only 
consider continuous wave experiments since this is what our instrument measures, where Eq. 
(1) is usually written with a jco/D(X) term on the right hand side, in our case we consider 
£0 = 0 giving us the form shown in Eq. (1). 

We employ the Born approximation by decomposing jtt°(r, A) and D°(r, A) , as the sum of a 
constant background absorption, jtt„(A), and a spatially varying perturbation Afi a (r,X) as well 
as constant background diffusion D(X) and a diffusion perturbation AD(r,X). To obtain a linear 
relationship between the measurements and the chromophore concentrations, we subtract Eq. 
(1) from the perturbed version of the diffusion equation which gives 

[V 2 +^(A)]4> i (r,A) = -A^(r,A)4>(r,A)-^ I yV.AD(r,A)V4> (3) 

where feg(A) = — v^ a (X)/D(X) and M 2 (r,A) = (v/D(A))A;U«(r,A). Assuming the availability 
of a Green's function, G(r,r',A) for the solution of Eq. (3) as is the case for an unbounded 
medium as well as range of bounded problems [26], we can change Eq. (3) into an integral 
Eq. [27] 

4>l(r ' A) = ~W) (/ A M(r',A)G(r,r',A)4>,(r',A) £ /r' 

+ J AD(r',A)VG(r f/ ,r„A)-V4> ; (r i ,r i .,Ayr') (4) 

where r c / is the location of the detector and (with a small abuse of notation) <I>, (r,r s ,A) is used 
here to denote the incident field at position r and wavelength A due to a delta-source located 
at r s . It should also be noted that we obtain this equation under the assumption that the total 
fluence rate, <E>, can be approximated as the incident field, <t>,, since <t> ( - 3> <i> s [2]. 

Equation (4) provides a linear relationship between the scattered fluence rate and the absorp- 
tion perturbation. To relate the scattered fluence rate to concentrations of chromophores, Aji a is 
decomposed as follows [28] 

N s 

AjU„(r,A) = £ e k (X)c k {r). (5) 

k=l 

where N s is the number of absorbing species for the problem under investigation, £a(A) is the 
extinction coefficient for the k' h species at wavelength A, and c^r) is the concentration of 
species k at location r. To obtain the fully discrete form of the Born model used in Section 4, 
we expand each Cfc(r) 

N, 

c*(r) = £ctjp(r) (6) 

7=1 

where c^j is the value of the concentration for species k in Vj, the j th "voxel". The (p(r) function 
is an indicator function where 

♦«-{« Vlv CO 

[0, lfr^V). 

After using Eqs. (5) and (6) in Eq. (4) and performing some algebra we obtain 
*v(A) = fl^y £ {G(r ch r j ,X)<P i {r j ,r S7 X)e i (X)c k j 

+ VG(r d ,ry,A)-V<J>Kr ; ,r,,A)AD,.(A)Y (8) 
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We approximate Eq. (4) as the value at the center of each pixel multiplied by the area or volume 
of each pixel or voxel, so in Eq. (8) we use a as the area of a pixel. This setup is illustrated in 
Fig. 1(a). 
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Fig. 1. (a) The setup of sources and detectors for simulation reconstructions. Same orienta- 
tion of axes is used for experimental data, (b) Definition of domains used for the parametric 
level-set methods. 

Setting up the linear algebraic structure associated with Eq. (8) we define Ck E M. Nv as the 
vector obtained by lexicographically ordering the unknown concentrations associated with the 
k lh chromophore and c^+i = A 1 ?' and <J>s(A) is the vector of observed scattered fluence rate 
associated with all source-detector pairs collecting data at wavelength A. Now, with Nx the 
number of wavelengths used in a given experiment, Eq. (8) is written in matrix-vector notation 
as 



> s (Ai)~ 




' ei(Ai)K? 


e 2 (Ai)K? . 
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Cl 
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<^ 4> v = Kc 



It should be noted in Eq. (9) 



trices K" and K^. The (m,j) th element of the K" matrix 
r ; -,A/)4>,-(r 7 -,r m 



that the matrix 
element of 
A/), where 1 



(9) 

has elements which are also the ma- 
is given by the prod- 
represents the m source-detector pair 
the matrix elements are given by 



Kf 



uct (v/£>(A,))G(r, 

and as before j represents the j"' voxel. For 
(v/D(A/)) VG(r m , ry, A) ■ V4>,- (ry , r m , A ) ADj (A ) . 

To evaluate the forward model for realistically sized problems, we compute the N\ matrices 
K" and K'' and store it along with the N% x N c extinction coefficients. This reduces the amount 
of memory required for the reconstruction. 



3. Parametric level-set method 

To counter the ill-posedness of the DOT problem we employ a Parametric Level-Set Method 
(PaLS). For the purpose of our paper we assume that all chromophore concentrations and dif- 
fusion coefficient perturbations are co-located. This choice is supported by reports in literature, 
where decrease in hemoglobin and water concentration along with scattering power are located 
at the cancer location and the lipid concentrations increase at the same location [29-31]. This 
means that the geometry of the anomaly in the medium is the same for all chromophores and 
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diffusion amplitude. The domain £2 C & represents the support of the objects of interest, and 
& represents the homogeneous background where the anomaly is located, shown in Fig. 1(b). 

Since the same support is used for each object of interest the characteristic function describ- 
ing the shape is defined as 

x^y)4 l if(x ' y)en ' do) 

Then each image to be reconstructed can be written as 

c k (x,y)=x(x,y)4+[i-x(x,y)H (ii) 

where k = 1,2, ...,N c + l. In this formulation the unknown values are the constant concentration 
values of the anomaly and background, c" k and c h k respectively. 

The characteristic function x(x,y) is defined to be the zero level set of a Lipschitz continuous 
object function G : & — ► K such that 0 > 0 in £l(x,y), 0 < 0 in and 0{x,y) = 0 in d£l. 
Using 0(x,y), x(x,y) is written as 

X{x,y)=H(ff(x,y)) (12) 

where H is the step function. In practice we use smooth approximations of the step function 
and the Dirac delta function denoted as H £ and 8 £ respectively where H £ is computed as 

if x > e, 

H e (x)={0 ifx<-e. 

'l + l + lMf)] if\x\<e. 

and 8 £ is computed as the derivative of H £ [32, 33]. As discussed in Section 1 we represent the 
object function 0(x,y) parametrically, so instead of using a dense collection of pixel or voxel 
values [34], we represent it by using basis functions 

L 

0(x 1 y) = '£a i p i (x,y) ( 13 > 

1=1 

where a,'s are the weight coefficients whereas pi(x,y) are the functions which belong to the 
basis set of = {pi,p2, ■ ■■,Pi}- Possible choices for the basis set include polynomial or 
radial basis functions. For the purpose of this paper we use Gaussian basis function. The width 
and number of the Gaussians determines how coarse or fine the reconstruction will be, where a 
choice of few basis functions will, on the one hand, result in a reduced number of unknowns, it 
will on the other hand, give a coarser estimation of the shape, which can be a problem for imag- 
ing finer more complex structures. In the DOT case, where the physics in the forward model 
will only allow for a coarse reconstruction of the underlying structure, the Gaussian approach 
is sufficient, especially for the relatively simple geometries and concentrations presented in 
this paper. When we will move on to a more complicated non-linear model, the choice of the 
number and type of basis functions will be more important and should be based on a rigid 
mathematical framework. This is discussed further in Section 8. 

In the PaLS formulation, all the parameters of the model are gathered in one vector 6 T = 
\c\ , , . . . , c( , a r ] where a = {a\, ...,ai] T . Now our forward model in Eq. (9) can be expressed as 

4>. S = K(0) =Kc(0) (14) 

The forward model has now been parametrized with a vector containing all of the unknowns, 
which are far fewer then what a pixel based method would attempt to estimate. 
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4. Inverse problem 



The inverse problem, that of using 4> s to recover the value of c, is solved as a Levenberg- 
Marquardt optimization problem of the form 

c = argmin||W(K(0)-4>. v )||^ (15) 

c 

The W matrix reflects the structure of the noise corrupting the data [2] . While a Poisson 
model is technically the most appropriate for DOT data, as is frequently done [35] we employ a 
Gaussian approximation in which independent, zero mean Gaussian noise is assumed to corrupt 
each datum. The reason for this is that with a sufficiently large number of detected photons, the 
Poisson statistics can be approximated by a Gaussian distribution [11]. Letting a} n denote the 
variance of the noise corrupting the m th elements of <t>, W is constructed as a diagonal matrix 
with 1/(7,,, the m th element along the diagonal. For the experimental and simulated data the 
variance is calculated from 

T , , SNR„, 

o£ = fi(m) 10 — nr. (16) 

where £2 (to) corresponds to the photon count for each source-detector pair. The SNR for each 
element of <t> is then calculated from 

SNR m = 101og 10 (£2(TO)/y£2R). (17) 

In experimental data y/Q,(m) is the standard deviation of the Poisson noise distribution. The 
minimization of the cost function is then achieved by the Levenberg-Marquardt algorithm. For 
that purpose an error vector is introduced 

e = W(K(0) -<!>). (18) 
which can be used to write the cost function in term of e as 



M(8) 



T 
£ £ 



(19) 



In order to employ the Levenberg-Marquardt algorithm, the calculation of the Jacobian ma- 
trix J is required. The Jacobian contains derivatives of e with respect to each element in the 
parameter vector 0 

de{8) 



d{cl,...c l b ,a}^ 

The solution is then obtained by updating 9 at each iteration as 6 
solution to the following linear system 

(J r J + pI)h = -fe with p > 0 



(20) 



0" + h where h is the 



(21) 



where I is the identity matrix, p is the damping parameter affecting the size and direction of 
h and found via and appropriate line search algorithm [36]. The stopping criteria used when 
iterating Eq. (21) is the discrepancy principle [37], in that the iterations are stopped when the 
norm of the residual has reached the noise level within a certain tolerance. 

When estimating the parametric vector, we employ a cyclic coordinate decent strategy [38] 
Essentially this is equivalent to estimating the shape only at even iterations and the concentra- 
tion values at odd iterations. This is repeated until stopping criteria is reached. This process 
is expressed in pseudo-code in Algorithm 1, where J v and J s denote the Jacobian strictly for 
the concentration value and shape, respectively, and T, represents a tolerance for the stopping 
critera. 
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Algorithm 1 Matlab-like code for estimating shape and concentration value 
while e < Ti do 
while e < T2 do 

(J v r J,. + pI)h rfl , !( « = -Jve 
end while 
while e < T3 do 

(JsJ s + pl)Khape = -jj£ 

end while 

9" +l = 9" + \hvalues'M s hape\ 

end while 



5. Simulation analysis 

The arrangement of sources and detectors with respect to the simulated turbid medium is dis- 
played in Fig. 1(a). This arrangement of sources and detectors is chosen to represent a common 
setup in optical mammography where the breast is compressed between two planes contain- 
ing sources and detectors [28]. The source-detector separations were set to 5 cm as is shown 
in Fig. 1(a). In simulations, we reconstruct concentration images of oxygenated and deoxy- 
genated hemoglobin, HbC>2 and HbR respectively, along with lipid and water concentration 
and diffusion amplitude. These chromophores are chosen since they mainly cause near-infrared 
absorption in breast tissue [39], and breast cancer tumours have been found to have higher 
HbO? and HbR concentrations than normal tissue [40]. 

In simulations, values for *P and b are obtained from [41] for the female breast. Values for jj. a 
are calculated from extinction coefficients, which are in the units cm^'/mM and are obtained 
from data tabulated by Scott Prahl [42]. The concentration in the simulated images are defined 
in units of millimolars or millimoles per liter, mM for HbC>2 and HbR. For water and lipid the 
concentrations are in percent by volume and the diffusion amplitude is measured in units of 
millimeter. The background has HbR concentration of 0.01 mM, HbC>2 concentration of 0.01 
mM, lipid concentration of 32%, water concentration of 13% and *P' is set to 1.6 mm. The 
target concentration of the object of interest is set to 0.015 mM, 0.012 mM, 50 %, 20 % and 
0.25 mm for Hb02, HbR, lipid, water and A 1 ?', respectively. 

The simulation set is created with all chromophore concentrations and diffusion perturbations 
in the same location with different target values. The co-located geometry is chosen since it is 
likely for chromophore concentrations to have similarly located geometries in the real breast 
[29, 30]. The ground truth images for simulations are shown in Fig. 6. Reconstruction is done 
for these images to explore effects of adding hyperspectral information to the problem, i.e. the 
improvement in quantitative accuracy and the reduction of crosstalk where a concentration of 
one chromophore creates a false concentration in an image for another chromophore as well 
as the performance of the shape based approach. The process is initialized with 21 Gaussian 
basis functions with width of approximately 8 pixels placed uniformly on a grid over the whole 
medium to be imaged. A representative image of the order of the basis functions is shown in 
Fig. 2(a). For all experiments presented in this paper, the a,'s weight coefficients are initialized 
to 1. 

To best understand the utility of a hyperspectral data set, we employ the Born model to 
generate simulated data. Though this may not be realistic, it allows us to avoid the confounding 
factor of model mismatch in evaluating the inversion method being considered in this paper. 
Moreover, the shortcomings of this approach are mitigated in Section 7.2, where we consider 
the processing of experimental data which, obviously, are not the product of the Born model. 
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Fig. 2. (a) Order of the basis functions used for the PaLS method, (b) Absorption spectra of 
the ink and dye solutions chromophores used in experimental measurements. Specifically 
chosen wavelengths are marked with an asterisk. 



In a bit more detail, the data we use for our simulation analysis are computed as 

* = Kc + n (22) 

where c are the simulated concentration images for all chromophores and diffusion amplitude, 
whereas n represents additive noise. Specifically, as in [2] n is a vector of zero mean, inde- 
pendent Gaussian random variables with variances a} n , defined in Eq. (16), chosen such that a 
pre-determined signal to noise ratio (SNR) is achieved. This SNR is calculated from Eq. (17). 

The reconstructed images are evaluated in three ways: through visual inspection, using mean 
square error (MSE) as a measure of overall quantitative accuracy for each chromophore, and 
examining the Dice coefficient to judge how well the concentrations are located. For the k th 
chromophore, the mean square error is computed by using the following Eq. 

MSE k = " , ,, 112 (23) 

Along with using MSE to judge the accuracy of the reconstruction, we also employ a Dice 
Coefficient to quantitatively analyze how well the shape based approach localizes the recon- 
struction [43,44]. If S is the reconstructed image and G is the ground truth created for each set, 
the Dice coefficient between S and G can be defined as 

\SllG\ contains all pixels that belong to the detected segment as well as the ground truth seg- 
ment, so that if S and G are equal the Dice coefficient is equal to one, indicating an accurate 
reconstruction. To compute the D(S,G) we use the characteristic function, %, which essentially 
works as a binary map of the reconstructed anomaly where the object of interest is represented 
by l's. 



6. Experimental analysis 

Physical measurements were performed in order to validate the simulation results, where the 
data used in this paper has previously been used in [10]. We use the same data again since we 
are using a completely new image formation process which results in a significant improvement 
in reconstructions. The background medium consists of water and milk in the ratio of 2:1, 
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respectively. Milk, with 2% fat, is used due to the similarities of the optical properties to human 
breast tissue. Black India ink and blue food dye were added to mimic tissue chromophores. The 
ink and dye are mixed into the background of milk and water to achieve jx a = 0.029 cm -1 , at 
600 nm, which is in the range of optical absorption of the female breast [45,46]. 

The absorption spectra for the ink and dye inclusions, shown in Fig. 2(b), have the most 
significant effect in the 450-700 nm range. These chromphores are chosen because the spectral 
shapes of their absorption are similar to those of Hb02 and HbR and have been widely used 
in literature [47,48]. Therefore, the 6 specifically chosen wavelengths were selected from this 
region where these chromophores had the most significant absorption. 

In order to obtain multi- and hyperspectral reconstruction values for jj, a and D(X) the back- 
ground has to be known. In the experimental measurements we assume constant scattering, 
therefore we are not trying to estimate the perturbation, AD(r,A), as was done in simulations. 
Therefore we have the unperturbed representaion of the reduced scattering coefficient, jit(, is 
given by 

*-!■(£) ■ (25) 

The diffusion coefficient relates to the reduced scattering coefficient through, D(X) = v/3fl' s . 
Phase, amplitude and average intensity data are obtained at two wavelengths using a frequency- 
domain tissue spectrometer to estimate the *P and b parameters in Eq. (25) as *F = 6.5 cm -1 
and b = 0.4. This allows us to extrapolate values for fi' s at any wavelength [49]. Determining 
spectrally extrapolated values for the absorption coefficient is harder. Since jj. a does not follow 
a law like /i s ', values are estimated using extinction coefficient data for ink, dye, milk and water. 
These extinction coefficient are measured in a standard spectrophotometer. 




Fig. 3. (a) Absorption spectra for the background, jl a , and the inclusion, il a + AjU a, in ex- 
perimental set 1, containing 10% ink and 90% dye. (b) Contrast between the background 
and the inclusion for experimental set 1 . 

In our experiments, two phantom inclusions, named set 1 and set 2, are created for different 
absorption contrasts relative to the background in the range of 3:1 to 1:1. The inclusion in set 
1 contains 10% ink and 90% dye and the inclusion for set 2 contains 70% dye and 30% ink. 
This contrast range is comparable to traditional tumor contrasts reported in literature, which 
have been close to 3:1 and lower [50]. Reconstructions are done for 126 wavelengths equally 
spaced over the whole spectrum and 6 specifically chosen wavelengths as A = [480, 550, 610, 
630, 650, 690] nm. The wavelengths are chosen around the isosbestic point, where the contrast 
between the chromophores is the highest and where each chromophore has highest absorption. 
The absorption spectra and the contrast over the spectrum for set 1 and set 2 are shown in Fig. 
3 and Fig. 4, respectively. 
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Wavelength [nm] 



560 580 600 620 640 660 
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Fig. 4. (a) Absorption spectra for the background, /l f( , and the inclusion, /J a + A/i fl , in 
experimental set 2, containing 70% ink and 30% dye. (b) Contrast between the background 
and the inclusion for experimental set 2. 

For both experimental sets 1 and 2, a single 25 cm long transparent cylindrical inclusion 
containing ink and dye solutions is placed in the background medium. Optical properties are 
assumed constant along the z-axis. For the light source, an arc lamp is used which delivers 
light with an average illumination power of 280 mW, that translates into a power density of 
3.96 W/cm 2 . For three different locations, at x ± 1 cm, a 5 mm diameter collection optical 
glass fiber bundle is placed on the opposite side of the inclusions (at a y-axis separation of 5 
cm). Experiments are made with the light source placed in succession at 8 positions with 1 cm 
increments for total of 24 source-detector pairs. To ensure that approximately the same amount 
of photons are collected for both hyperspectral and multispectral reconstructions, two exposure 
times are used for the CCD camera, a longer one of 10 s for the 6 wavelength case and 500 ms 
for the 126 wavelength case. Since the goal of this paper is to demonstrate the improvement 
of including hyperspectral information, we present an ideal case where the signal to noise ratio 
is large, thereby providing a best-case scenario for the few wavelength reconstruction against 
which we compare our approach as well as using realistic absorption contrasts for the inclu- 
sions. Further details on the experimental setup can be found at [10]. 

In our experimental setup, we measure the incident field before the perturbation is put into 
the medium. Then the scattered field is computed as a dataset that has the unperturbed dataset 
subtracted from it. For in vivo measurements it is possible to use a priori structural information 
from other modalities, e.g. MRI, to estimate the incident field by determining the optical prop- 
erties of the assumed piecewise constant chromophore distribution over these segments [51]. 

A comparison of the absolute concentrations, c} and relative concentration, c- to target con- 
centration values is done to test the accuracy of the reconstructions. The relative concentrations 
for ink are calculated as 

C l,k =^ink/ i^ink + Cctye) (26) 

and similarly for dye [52,53]. The relative concentration is calculated from the peak concentra- 
tion value in each reconstruction. This allows us to inspect how well our approach manages to 
separate and estimate each species of chromophores in the process. 

7. Results 

7.1. Simulations 

In Fig. 5 reconstruction results using the pixel based method are shown for 8 wavelengths, A = 
[660, 734, 760, 808, 826, 850, 930, 980] nm and hyperspectral reconstruction using 176 wave- 
lengths, which are equally spaced over the 650-1000 nm range. In the 8 wavelength case the 
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2468 10 2468 10 2468 10 



Fig. 5. Reconstruction using a pixel based method. Middle column of images are generated 
with 8 wavelengths and rightmost images are generated with 176 wavelengths. From top 
to bottom the rows show HbC>2, HbR, lipid, water and diffusion amplitude, respectively. 
Concentration units are in mM. 

Table 1. The MSE is compared for each chromphore for multiple wavelength choices. In 
each case the reconstructions are done with equally spaced wavelengths over the spectrum 
except for the 8 wavelength case. 



Pixel based method 


#A 


MSE Hb0 2 


MSE HbR 


MSE Lipid 


MSE H 2 0 


MSE D 


8 


0.075 


0.030 


0.048 


0.010 


0.052 


176 


0.062 


0.021 


0.034 


0.015 


0.030 


PaLS method 


# X 


MSE Hb0 2 


MSE HbR 


MSE Lipid 


MSE H 2 0 


MSE D 


8 


0.07 


0.03 


0.12 


0.06 


0.08 


176 


0.02 


0.008 


0.01 


0.02 


0.01 



6 first wavelengths are optimally chosen according to [54] with two wavelengths added where 
water and lipids have peak absorption. Reconstructed images created with the PaLS method are 
shown in Fig. 6. In simulations the SNR is set to 30 dB, as it is defined by Eq. (16) and Eq. (17). 
When comparing the pixel based reconstruction in Fig. 5 to the PaLS reconstruction in Fig. 6, it 
is evident that the PaLS method provides superior reconstructions. Examining the PaLS results, 
the 8 wavelength case shows reasonable accuracy along the x axis but rather diffuse results in 
y. We also see noticeable artifacts in the reconstructions. Considering the concentration values, 
the values for HbC<2, HbR and water concentration come close to the actual value. Moving 
to hyperspectral information, the reconstruction becomes more accurate, estimating the shape 
close to the ground truth. It should also be noted that the runtime for each reconstruction for 
the PaLS method is significantly shorter compared to the pixel-based method. A PaLS recon- 
struction takes around 30 seconds, which is 3-4 times faster than a pixel-based method. Addi- 
tionally, we do not employ any regularization parameters, freeing us from finding the optimal 
reconstruction using regularization. This is a major improvement in moving from a pixel-based 
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Table 2. D(S,G) is compared for each chromphore for multiple wavelength choices. In 
each case the reconstructions are done with equally spaced wavelengths over the spectrum 
except for the 8 wavelength case. D(S, G) is calculated comparing 80% of the target peak 
to the reconstruction. 



Pixel based method 


# X 


D(S,G) Hb0 2 


D(S,G) HbR 


D(S,G) Lipid 


D(S,G) H 2 0 


D(S,G) D 


8 


0.12 


0.088 


0.089 


0.65 


0.8 


176 


0.554 


0.1085 


0.043 


0.41 


0.09 




PaLS method 




#A 


D{S,G) 


8 


0.60 


126 


0.99 




'\ •HQ- '\ ID :;! ■] P 

2 4 6 8 10 2468 10 2468 10 



Fig. 6. Reconstruction using the PaLS method. Middle column of images are generated 
with 8 wavelengths and rightmost images are generated with 176 wavelengths. From top 
to bottom the rows show HbC>2, HbR, lipid, water and diffusion amplitude, respectively. 
Concentration units are in mM. 



approach to the PaLS method. 

The comparison of the Dice coefficient between the PaLS method and pixel-based is tricky, 
since for the pixel-based method the Dice coefficient is plotted as a function of a threshold. 
This threshold is required to create a binary map of the location on the anomaly. If the threshold 
is chosen to only leave extreme peak concentration values in each image, the Dice coefficient 
would be low due to edge artifacts as in Fig. 7(b). Therefore, in simulations we compare D(S, G) 
for the pixel based reconstructions using a threshold of 80% to D(S, G) of the PaLS reconstruc- 
tions. The improvement of the PaLS method is confirmed quantitatively through D(S,G) and 
MSE displayed in Table 2 and Table 1, respectively. The Dice coefficient, shown in Table 2, 
gives a clear view of how the shape estimation improves by added wavelengths, where D(S, G) 
approaches 1 for the hyperspectral case and the PaLS method shows superior performance in 
the MSE values. 
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7.2. Experimental validation 




cm 



(a) 10% ink and 90% dye, 6 wavelengths used. (b) 70% ink and 30% dye, 6 wavelengths used. 




(c) 10% ink and 90% dye, 126 wavelengths used. (d) 70% ink and 30% dye, 126 wavelengths used. 

Fig. 7. Pixel reconstruction from both experimental sets, set 1 containing 10% ink and 90% 
dye and set 2 70% ink and 30% dye. 

Pixel based reconstructions of absolute concentrations for both experimental sets are shown 
in Fig. 7 and PalS approach reconstructions in Fig. 8. As expected, the hyperspectral informa- 
tion provides improved reconstruction for both the pixel-based and PaLS methods. The pixel 
based results had previously been reported in [10]. Focusing on the PaLS methods, it is evident 
that forming the reconstruction with shape-based constraints yields improved results. The esti- 
mation of relative concentrations and MSE of the absolute values are examined in Table 4 and 
Table 5 for the pixel-based and PaLS method, respectively. The relative concentration values are 
better estimated in both cases using the PaLS method, although the hyperspectral method does 
not show significant improvement for experimental set 2, which was also the case for the pixel 
based method. Examining the images along with the MSE values for experimental set 1, Fig. 
7-8(a) and (c), it is noticeable how the reconstruction does not resolve the structure particurarly 
well along the x axis. This is somewhat unexpected since in DOT resolving depth information, 
on the y axis, is usually the more difficult problem. This is noticeable for both the pixel based 
and PaLS methods, although the PaLS method outperforms the pixel based method, especially 
in removing edge artifacts. This streaking in the x direction is most likely a combination of how 
the Gaussian basis are placed within the imaging medium, and measurement error in placing 
the source and detectors when taking the reference measurement. 

For both experimental sets, the PaLS method resolves the location and the shape of the in- 
clusion more accurately, which is verified by the calculation of the Dice coefficient shown in 
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Table 3, the improvement is notable when compared to the pixel-based reconstruction. As dis- 
cussed in Section 7.1a choice of a threshold is needed to compare D(S, G) between pixel based 
reconstructions and the PaLS method. For the experimental reconstructions we use a threshold 
of 50% to compare D(S,G) of the PaLS reconstructions. This demonstrates the usefulness of 
the PaLS method correctly and accurately localizing the anomaly. The PaLS method does very 
well with eliminating edge artifacts that were severe when doing pixel-based reconstructions 
for the same data set. These effects are very noticeable in Fig. 7(b) and (d), where, especially 
in the multispectral case, the edge artifacts were significant. Comparing that to the same data 
in Fig. 8(b) and (d) it is obvious that the improvement is significant. 



Table 3. D(S,G) is compared for each chromphore for multiple wavelength choices. In 
each case the reconstructions are done with equally spaced wavelengths over the spectrum 
except for the 6 wavelength case where we use optimally chosen wavelengths. D(5, G) is 
calculated comparing the half maximum of the target peak to the reconstruction. 



Pixel based method 




D{S,G) Set 1 


D(S,G) Set 2 


# A 


Ink 


Dye 


Ink 


Dye 


6 


0.143 


0.113 


0.139 


0.145 


126 


0.142 


0.114 


0.145 


0.140 



PaLS method 


# A 


D(S,G) Set 1 


D(S,G) Set 2 


6 


0.27 


0.33 


126 


0.37 


0.80 



Table 4. Comparison of c, and c[ to target concentration values for experimental results, for 
the pixel-based method. Best performance is highlighted in bold. 



Experimental set 1, 10% ink and 90% dye 


Fig. 


#A 


Species 


Cj [%] 


cj [%] 


MSE 


7(a) 


6 


Ink 


1 


4 


1.8 


7(a) 


6 


Dye 


27 


96 


1.3 


7(c) 


126 


Ink 


17 


16 


2.8 


7(c) 


126 


Dye 


88 


84 


1.2 



Experimental set 2, 70% ink and 30% dye 


Fig. 


#A 


Species 


c, [%] 


<? [%i 


MSE 


7(b) 


6 


Ink 


56 


82 


1.8 


7(b) 


6 


Dye 


12 


18 


1.0 


7(d) 


126 


Ink 


65 


61 


1.4 


7(d) 


126 


Dye 


41 


39 


2.0 



Table 5. Comparison of c, and cj to target concentration values for experimental results, for 
the PaLS method. Best performance is highlighted in bold. 



Experimental set 1, 10% ink and 90% dye 


Fig. 


#A 


Species 


c, [%] 


< [%] 


MSE 


8(a) 


6 


Ink 


4.8 


21.0 


1.2 


8(a) 


6 


Dye 


17.9 


79.0 


0.9 


8(c) 


126 


Ink 


5.8 


7.7 


1.1 


8(c) 


126 


Dye 


69.0 


92.3 


0.8 



Experimental set 2, 70% ink and 30% dye 


Fig. 


#A 


Species 


c, [%] 


cf [%] 


MSE 


8(b) 


6 


Ink 


38.3 


80.0 


1.1 


8(b) 


6 


Dye 


9.6 


20.0 


0.8 


8(d) 


126 


Ink 


27.6 


81.0 


0.6 


8(d) 


126 


Dye 


6.4 


19.0 


0.7 



8. Conclusion 

In this paper, using simulations and experimental measurements we have shown that the PaLS 
method provides more accurate estimation of chromophore concentrations than a regularized 
pixel-based inversion scheme. Hyperspectral information results in improved performance in 
terms of both MSE and spatial localization as measured using the Dice coefficient. The para- 
metric approach is shown to give significant improvements to image reconstruction decreasing 
run time of the iterative process and increasing the quality of reconstructed images. The PaLS 
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cm 



(a) 10% ink and 90% dye, 6 wavelengths used. (b) 70% ink and 30% dye, 6 wavelengths used. 




cm 



Reconstructed image of Dye 



Reconstructed image of Dye 



(c) 10% ink and 90% dye, 126 wavelengths used. (d) 70% ink and 30% dye, 126 wavelengths used. 

Fig. 8. PaLS Reconstruction from both experimental sets, set 1 containing 10% ink and 
90% dye and set 2 70% ink and 30% dye. 



method is also easily expandable to more complicated problems where multiple geometries 
need to be considered. 

Physical measurements were also performed to demonstrate these advantages for actual 
measurement data. Although exact concentration values were not achieved, there is a notable 
improvement associated with hyperspectral information in conjunction with the PaLS method. 
Additionally, improved localization of inclusions was observed for both sets when using hyper- 
spectral information. This emphasizes the advantage of hyperspectral information when doing 
reconstructions for more than one chromophore. 

Based on the results in this paper, we will be extending the work to address more realistic, 
clinical conditions. We will first implement a fully 3D model employing boundary conditions 
and consider an inhomogeneous background more similar to what is found in breast imaging. In 
that setting, estimating multiple level sets for different geometries could prove useful, especially 
to estimate different regions of the heterogeneous background such as adipose and fibroglan- 
dular tissue. As mentioned in Section 3, it is also important to generate a rigid framework to 
choose different types, sizes and shapes of basis functions to generate the best reconstruction. 
To be able to estimate all shapes possible, we aim to increase the number of basis functions 
to include different types. To avoid over complicating the image reconstruction with a high 
number of basis functions we aim to pose the image reconstruction in a compressed sensing 
framework where few optimal basis functions estimate a complex shape. Furthermore, relating 
to the development of a framework for different basis functions, we will develop the Leven- 
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berg Marquardt algorithm to have optimally chosen stopping criteria and step size changes. 
This should increase robustness of our method and ensure accuracy of estimation for different 
situations and settings in diffuse optical tomography. 
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